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Abstract 

A transonic unsteady aerodynamic and aeroelasticity code called CAP-TSD has been developed 
for application to realistic aircraft configurations. The code permits the calculation of steady 
and unsteady flows about complete aircraft configurations for aeroclastic analysis in the flutter 
critical transonic speed range. The CAP-TSD code uses a time-accurate approximate 
factorization (AF) algorithm for solution of the unsteady transonic small-disturbance potential 
equation. The paper gives an overview of the CAP-TSD code development effort and presents 
results which demonstrate various capabilities of the code. Calculations are presented for 
several configurations including the General Dynamics one-ninth scale F-16C aircraft model 
and the ONERA M6 wing. Calculations are also presented from a flutter analysis of a 45° 
sweptback wing which agree well with the experimental data. The paper presents descriptions 
of the CAP-TSD code and algorithm details along with results and comparisons which 
demonstrate these recent developments in transonic computational aeroelasticity. 

Notation 

b 0 reference length, cr/2 

c, cr airfoil chord and wing reference chord, respectively 

Cp pressure coefficient 
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k reduced frequency, tocr/21) 

M freestream Mach number 

t time, nondimensionalized by freestream speed and wing reference chord 

U freestream speed 

ao, ai mean angle of attack and amplitude of pitch oscillation, respectively 
y ratio of specific heats 

At nondimensional time step 

r\ fractional semispan 

p ratio of wing mass to mass of air in the truncated cone that encloses the wing 

p freestream flow density 

<t> disturbance velocity potential 

o), cou angular frequency and natural frequency of the first torsion mode, respectively 

Subscripts 
t tail 

w wing 


Introduction 

Presently, considerable research is being conducted to develop finite-difference computer 
codes for calculating transonic unsteady aerodynamics for aeroelastic applications [1]. These 
computer codes are being developed to provide accurate methods of calculating unsteady airloads 
for the prediction of aeroelastic phenomena such as flutter and divergence. For example, the 
CAP-TSD [2] unsteady transonic small-disturbance (TSD) code was recently developed for 
transonic aeroelastic analyses of complete aircraft configurations. The name CAP-TSD is an 
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acronym for Computational Aeroelasticity Program - Transonic Small Disturbance. The new 
code permits the calculation of unsteady flows about complete aircraft for aeroelastic analysis 
in the flutter critical transonic speed range. The code can treat configurations with arbitrary 
combinations of lifting surfaces and bodies including canard, wing, tail, control surfaces, tip 
launchers, pylons, fuselage, stores, and nacelles. Steady and unsteady pressure comparisons 
were presented for numerous cases which demonstrated the geometrical applicability of CAP- 
TSD [2-3]. These calculated results were generally in good agreement with available 
experimental pressure data which validated CAP-TSD for multiple component applications with 
mutual aerodynamic interference effects. Aeroelastic applications of CAP-TSD were presented 
by Cunningham, et al. [4] and Bennett, et al. [5] for simple well-defined wing cases. The cases 
were selected as a first step toward performing aeroelastic analyses for complete aircraft 
configurations. The calculated flutter boundaries compared well with the experimental data for 
subsonic, transonic, and supersonic freestream Mach numbers, which gives confidence in CAP- 
TSD for aeroelastic prediction. 

The CAP-TSD code uses a time-accurate approximate factorization (AF) algorithm recently 
developed by Batina [6] for solution of the unsteady TSD equation. The AF algorithm involves a 
Newton linearization procedure coupled with an internal iteration technique. The algorithm was 
shown to be efficient for application to steady or unsteady transonic flow problems. It can 
provide accurate solutions in only several hundred time steps, yielding a significant 
computational cost savings when compared to alternative methods. For reasons of practicality 
and affordability, an efficient algorithm and a fast computer code are requirements for realistic 
aircraft applications. 

Recently, several algorithm modifications have been made which have improved the stability 
of the AF algorithm and the accuracy of the results [7,8], These algorithm modifications 



4 


include: (1) an Engquist-Osher (E-0) type-dependent switch to more accurately and 

efficiently treat regions of supersonic flow, (2) extension of the E-0 switch for second-order- 
accurate spatial differencing in supersonic regions to improve the accuracy of the results, (3) 
nonisentropic effects to more accurately treat cases with strong shocks, and (4) nonreflecting 

i far field boundary conditions for more accurate unsteady applications. The work has been a 

major research activity over the past two and one-half years within the Unsteady Aerodynamics 
Branch at NASA Langley Research Center. The purpose of the paper is to give an overview of the 
CAP-TSD code development effort and report on the recent algorithm changes and code 
improvements. The paper documents these developments and presents results which 
demonstrate these recent advances in transonic computational aeroelasticity. 

Transonic Small Disturbance Equation 

| The flow is assumed to be governed by the general frequency modified TSD potential equation 

which may be written as 

M 2 », + 2* x ), = [(1 - M\ + F<f 2 + + (<(. y + + (<y x ( , , 

Several choices are available for the coefficients F, G, and H depending upon the assumptions 
used in deriving the TSD equation. For transonic applications, the coefficients are herein 
defined as 

i 

F = - — (Y + 1) M 2 » G= — (y-3)M 2 , H = -(y-l)M" (2) 

| 

The linear potential equation is recovered by simply setting F, G, and H equal to zero. 
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Entropy and vorticity modifications to TSD theory, to treat cases with strong shock waves, 
have been developed as described by Batina [8]. These modifications include: (1) an alternative 
streamwise flux in the TSD equation which was derived by an asymptotic expansion of the Euler 
equations, (2) a modified velocity vector defined as the sum of potential and rotational parts 
which in turn modified the streamwise flux, and (3) the calculation and convection of entropy 
throughout the flow field. The modified theory includes the effects of entropy and vorticity 
while retaining the relative simplicity and cost efficiency of the TSD formulation [8]. 

Approximate Factorization Algorithm 

A time-accurate approximate factorization algorithm was developed [6-8] to solve the 
unsteady TSD equation including entropy and vorticity effects. In this section, the AF algorithm 
is briefly described. 

General Description 

The AF algorithm consists of a Newton linearization procedure coupled with an internal 
iteration technique. For unsteady flow calculations, the solution procedure involves two steps. 
First, a time linearization step is performed to determine an estimate of the potential field. 
Second, internal iterations are performed to minimize linearization and factorization errors. 
Specifically, the TSD equation is written in general form as 

R(<t>" +1 ) = 0 (3) 

where <j> n+1 represents the unknown potential field at time level (n+1). The solution to Eq. (3) 
is then given by the Newton linearization of Eq. (3) about <>' 
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R (0*) + (—).* A<t> = 0 (4) 

In Eq. (4), 0* is the currently available value of 0 n+1 and A0 = 0 n+1 - 0*. During convergence 
of the iteration procedure, A0 will approach zero so that the solution will be given by 0 n+1 = 
0*. In general, only one or two iterations are required to achieve acceptable convergence. 

Mathematical Formulation 

The AF algorithm is formulated by first approximating the time-derivative terms by 
second-order-accurate finite-difference formulae. The TSD equation is rewritten by 

substituting 0 = 0* + A0 and neglecting squares of derivatives of A0 (which is equivalent to 
applying Eq. (4) term by term). The resulting equation is then rearranged and the left-hand 
side is approximately factored into a triple product of operators yielding 

L L L A0 = - a R(0*, <J>", 0"' 1 , <|> n ' 2 ) ( 5 ) 

where the operators L^, Lq, U; and residual R were derived and presented in Ref. [6]. In Eq. 

(5) ct is a relaxation parameter which is normally set equal to 1.0. To accelerate convergence 
to steady-state, the residual R may be over-relaxed using a > 1. Equation (5) is solved using 

three sweeps through the grid by sequentially applying the operators L4, U|, L; as 


£ - sweep: A0 = - a R 


(6a) 
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T| - sweep: 

L A<)> = A<t> 
n 

(6b) 

C - sweep: 

L A0 = A<j> 

(6c) 


The AF algorithm uses an Engquist-Osher (E-O) type-dependent switch to change the spatial 
differencing from central differencing in regions of subsonic flow to upwind differencing in 
regions of supersonic flow. This, of course, allows for the correct numerical description of the 
physical domain of dependence. The E-0 switch is based on sonic reference conditions and does 
not admit expansion shocks as part of the solution. Use of the E-0 switch also generally 
increases computational efficiency because of the larger time steps which may be taken. 
Furthermore, the E-0 switch of the AF algorithm has been recently extended for second-order 
spatial accuracy in supersonic regions of the flow. Details of these developments are reported 
by Batina [7). 

Boundary Conditions 

The flow tangency boundary conditions are imposed along the mean plane of the respective 
lifting surfaces and the wakes are assumed to be planar extensions from the trailing edges to the 
downstream boundary of the finite-difference grid. The numerical implementation of these 
conditions [2] allows for coplanar as well as non-coplanar combinations of horizontal (canard, 
wing, horizontal tail, launchers) and vertical (pylons, vertical tail) surfaces. Bodies such as 
the fuselage, stores, and nacelles are treated using simplified boundary conditions on a 
prismatic surface rather than on the true surface [2]. The method is consistent with the small- 
disturbance approximation and treats bodies with sufficient accuracy to obtain the correct 
global effect on the flow field without the use of special grids or complicated coordinate 



8 


transformations. This type of modeling is similar to that of Boppe and Stern [9] where good 
agreement was shown in comparison with experimental data for configurations with a fuselage 
and flow-through nacelles. 

For unsteady applications, nonreflecting far field boundary conditions similar to those of 
Whitlow [10] are used. The nonreflecting conditions absorb most of the waves that are incident 
on the grid boundaries and consequently allow the use of smaller computational grids. These 
boundary conditions are consistent with the AF solution procedure and are described in more 
detail by Batina [7]. 


Aeroelastic Solution 

In this section the aeroelastic computational procedures are described including the 
equations of motion and the time-marching solution. 

Eflya t jQ_n_s._Q f MotiQn 

The aeroelastic equations of motion are based on a right-hand orthogonal coordinate system 
with the x-direction defined as positive downstream and the z-direction positive upward [4]. 
The presentation herein is limited to the case of an isolated wing with motion in the z-direction 
from an undisturbed position in the z = 0 plane. The equations of motion may be written as 

Mq + Cq + Kq = Q ( 7 j 

where q is a vector of generalized displacements, M is the generalized mass matrix, C is the 
damping matrix, K is the stiffness matrix, and Q is the vector of generalized aerodynamic 


forces. 
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Time-Marchina Solution 

The aeroelastic solution procedure for integrating Eq. (7) is similar to that described by 
Edwards, et al. [11]. Reference [11] describes for a two-dimensional, two-degree-of-freedom 
system an aeroelastic solution in terms of a state equation formulation. By a parallel 
formulation, a linear state equation is developed from Eq. (7) which is solved numerically using 
the modified state-transition matrix integrator of Ref. [11]. This integrator was shown to be 
superior to six alternative structural integration algorithms [12]. 

For aeroelastic analysis, two steps are generally required in performing the calculations. In 
the first step, the steady-state flow field is calculated to account for wing thickness, camber, 
and mean angle of attack thus providing the starting flow field for the aeroelastic analysis. The 
second step is to prescribe an initial disturbance to begin the structural integration. 
Disturbance velocities in one or more modes, rather than displacements, have been found to be 
distinctly superior in avoiding nonphysical, strictly numerical transients and their possible 
associated instabilities. In determining a flutter point, the freestream Mach number M and the 
associated freestream speed U are usually held fixed. A judiciously chosen value of the dynamic 

pressure pU 2 /2 is used to compute the free decay transients. These resulting transients of the 
generalized coordinates are analyzed for their content of damped or growing sine-waves, with 
the rates of growth or decay indicating whether the dynamic pressure is above or below the 
flutter value. This analysis then indicates whether to increase or decrease the value of dynamic 
pressure in subsequent runs to determine a neutrally stable result. 

CAP zI SP .Code 

The AF algorithm is the basis of the CAP-TSD code for transonic unsteady aerodynamic and 
aeroelastic analysis of realistic aircraft configurations. The present capability has the option of 


half-span modeling for symmetric cases or full-span modeling to allow the treatment of 
antisymmetric mode shapes, fuselage yaw, or unsymmetric configurations such as an oblique 
wing or unsymmetric wing stores. In the present coding of the AF algorithm, the time 
derivatives are implemented for variable time stepping to allow for step-size cycling to 

accelerate convergence to steady state. Also, since the L^, L^, and operators only contain 

derivatives in their respective coordinate directions, all three sweeps of the solution procedure 
have been fully vectorized. 


Results and Discussion 

Results are presented for several configurations including the General Dynamics one-ninth 
scale F-16C aircraft model [13,14], the ONERA M6 wing [15], and a 45° sweptback wing [16, 
17], to demonstrate various capabilities of the CAP-TSD code. 

General Dynamics F- 16 C Aircraft Model Results 

To demonstrate application of CAP-TSD to a realistic aircraft configuration, results were 
obtained for the General Dynamics one-ninth scale F-16C aircraft model. Shown in Fig. 1 are 
the F-16C components that are modeled using CAP-TSD. The F-16C is modeled using four 
lifting surfaces and two bodies. The lifting surfaces include: (1) the wing with leading and 
trailing edge control surfaces, (2) the launcher, (3) a highly-swept strake, aft strake, and 
shelf surface, and (4) the horizontal tail. The bodies include: (1) the tip missile, and (2) the 
fuselage. In these calculations, the freestream Mach number was M = 0.9 and the aircraft was 
at 2.38° angle of attack. Also, the leading edge control surface of the wing was deflected upwards 
2° for comparison with the experimental steady pressure data of Ref. [14]. These steady 
pressure comparisons are made to assess the accuracy of CAP-TSD for complete airplane 


applications. There are no unsteady experimental data to validate the CAP-TSD code for time 
accurate F-16C calculations. Nonetheless, an unsteady calculation was performed to 
demonstrate the time-accurate capability. For simplicity, the results were obtained for a rigid 
pitching motion where the entire F-16C aircraft was forced to oscillate about the model moment 
reference axis at a reduced frequency of k = 0.1. The oscillation amplitude was chosen as ai = 
0.5°, and 300 steps per cycle of motion were computed corresponding to At = 0.1047. Parallel 
results were also obtained for the wing alone to investigate the aerodynamic interference effects 
of the additional aircraft components on wing unsteady pressures. 

Steady pressure comparisons for the F-16C aircraft model are presented in Fig. 2 for three 
span stations of the wing and one span station of the tail. For this case, there is a moderately 
strong shock wave on the upper surface of the wing and the CAP-TSD pressures agree well with 
the experimental pressures. For the tail, the flow is predominantly subcritical and the CAP- 
TSD pressures again agree well with the experimental data. 

Unsteady pressure distributions along the wing and tail upper surfaces are shown in Fig. 3 
for the entire F-16C aircraft undergoing a rigid pitching motion. These unsteady pressure 
results are presented at the same span stations as the steady-state results (Fig. 2). Two sets of 
calculated pressures are compared corresponding to complete airplane and wing alone modeling. 
As shown in Fig. 3, there is a relatively large shock pulse in the real part of the wing upper 
surface pressures due to the motion of the shock wave. This shock pulse is of larger magnitude 
and is located further downstream in the complete airplane model. These features are attributed 
to a stronger steady-state shock on the upper surface of the wing produced by the accelerated 
flow about the fuselage and the launcher/tip missile. For the tail, the unsteady pressures are 
relatively small in comparison with the wing pressures and thus were plotted on an expanded 
scale. The tail is located considerably aft of the pitch axis and thus its motion is plunge 


1 2 

dominated which results in smaller airloads for the low value of k considered. Furthermore, 
these pressures are nearly 90° out of phase with the aircraft motion since the real components 
are small compared to the imaginary components. 

ONERA M6 Wing Results 

To test the entropy and vorticity modifications to TSD theory, applications were performed 
for the ONERA M6 wing [15]. The M6 wing has an aspect ratio of 3.8, a leading edge sweep angle 
of 30°, and a taper ratio of 0.562. The airfoil section of the wing is the ONERA "D" airfoil 
which is a 10% maximum thickness-to-chord ratio symmetric section. Pressures were 
calculated at M = 0.92 with the wing at 0° mean angle of attack. These conditions correspond to 
an AGARD test case for assessment of inviscid flow field methods [15] and were selected for 
comparison with the tabulated Euler results of Rizzi contained therein. 

Calculations were performed using: (a) unmodified TSD theory and (b) TSD with entropy 
and vorticity effects. Steady pressure distributions along three span stations (fj = 0.08, 0.47, 
and 0.82) of the wing are presented in Fig. 4 from both solutions. For this case, the flow is 
symmetric about the wing with shocks on the upper and lower surfaces. As shown in Fig. 4(a), 
the results from the unmodified TSD theory compare well with the Euler results in predicting 
the leading edge suction peak and the overall pressure levels. However, the shock is located too 
far aft and is too strong outboard near the tip in comparison with the Euler calculation. When 
the entropy and vorticity modifications are included in the calculation, the shock is displaced 
forward from the previous solution, as shown in Fig. 4(b). Here the shock location and shock 
strength are in very good agreement with the Euler results at all three span stations on the 
wing. Consequently, the steady pressure distributions from the modified TSD theory now 
compare very well with the Euler pressures. 


Wing Flutter Results 


To assess the CAP-TSD code for flutter applications, a simple well-defined wing case was 
selected as a first step toward performing aeroelastic analyses for complete aircraft 
configurations [4,5]. The wing being analyzed is a semispan wind-tunnel-wall-mounted model 
that has a quarter-chord sweep angle of 45°, a panel aspect ratio of 1 .65, and a taper ratio of 
0.66 [16], The wing is an AGARD standard aeroelastic configuration which was tested in the 
Transonic Dynamics Tunnel (TDT) at NASA Langley Research Center [17], A planview of the 
wing is shown in Fig. 5. The wing has a NACA 65A004 airfoil section and was constructed of 
laminated mahogany. In order to obtain flutter for a wide range of Mach number and density 
conditions in the TDT, holes were drilled through the wing to reduce its stiffness. To maintain 
the aerodynamic shape of the wing, the holes were filled with a rigid foam plastic. The wing is 
modeled structurally using the first four natural vibration modes which are illustrated in Fig. 
6. These modes which are numbered 1 through 4 represent first bending, first torsion, second 
bending, and second torsion, respectively, as determined by a finite element analysis. The 
modes have natural frequencies which range from 9.6 Hz for the first bending mode to 91.54 Hz 
for the second torsion mode. 

Flutter calculations were performed for the 45° sweptback wing using CAP-TSD to assess 
the code for aeroelastic applications. Two sets of results are presented corresponding to: (1) 
using the linear potential equation (F = G = H = 0) and modeling the wing aerodynamically as a 
flat plate (zero thickness) and (2) using the complete (nonlinear) TSD equation and including 
wing thickness. The first set of results allows for direct comparison with parallel linear 
theory calculations performed using the FAST subsonic kernel function program [18]. The 
second set of results more accurately models the wing geometry as well as the flow physics. All 
of the results are compared with the experimental flutter data of Yates, et al. [17] which spans 
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the range, 0.338 < M < 1 .141 . 

Comparisons of flutter characteristics from the linear calculations with the experimental 
data are given in Fig. 7. The flutter speed index (defined as U/(bo ua Vu)) as a function of 

freestream Mach number is shown in the figure. The experimental flutter data defines a typical 
transonic flutter "dip” with the bottom near M = 1.0 for this case. The bottom of the dip in 
flutter speed index was defined by the approach to the M = 1.072 flutter point during the wind 
tunnel operation. Results from the CAP-TSD (linear) code are presented at twelve values of M 
covering the entire Mach number range over which the flutter data was measured. Results from 
the FAST program are presented for the limited range 0.338 < M ^ 0.96 since the method is 
restricted to subsonic freestreams. Overall, the linear CAP-TSD results compare well with the 
experimental data for subsonic as well as supersonic Mach numbers. Note that the subsonic 
FAST results are also in good agreement with the data. Such a result is not unexpected for this 
very thin wing of moderate sweep and taper at zero angle of attack. It does indicate that the wing 
properties are well-defined for benchmark purposes. 

In the subsonic Mach number range, the CAP-TSD and FAST calculations predict a slightly 
unconservative flutter speed, except at M = 0.338, by as much as 2% (Fig. 7) in comparison 
with the experimental data. In general though, the linear CAP-TSD results agree well with the 
FAST results. The good agreement in this three-way correlation between experiment, linear 
theory, and CFD flutter results gives confidence in the CAP-TSD code for flutter prediction. 

Comparisons of flutter characteristics from the linear and nonlinear CAP-TSD calculations 
with the experimental data are also given in Fig. 7. Three flutter points are plotted from the 
nonlinear CAP-TSD calculations corresponding to M = 0.678, 0.901, and 0.96. Comparisons 
between the two sets of CAP-TSD results show differences due to wing thickness and nonlinear 
effects. With increasing Mach number these differences become larger. For example, at M = 



0.678, 0.901, and 0.96, the flutter speed index decreased by 1%, 5%, and 19%, respectively. 
The decrease in flutter speed at M = 0.901 is largely due to including wing thickness since there 
are no supersonic points in the flow at this condition. The decrease in flutter speed at M = 0.96 
is attributed to both wing thickness and nonlinear effects since an embedded supersonic region of 
moderate size was detected in the wing tip region. The nonlinear CAP-TSD results at both M = 
0.901 and 0.96 are slightly conservative in comparison with the experimental flutter speed 
index values. Nonetheless, the nonlinear CAP-TSD flutter results compare favorably with the 
experimental data, which is the first step toward validating the code for general transonic 
aeroelastic applications. 


Concluding Remarks 

A transonic unsteady aerodynamic and aeroelasticity code called CAP-TSD has been developed 
for application to realistic aircraft configurations. The code permits the calculation of unsteady 
flows about complete aircraft configurations for aeroelastic analysis in the flutter critical 
transonic speed range. The CAP-TSD code uses a time-accurate approximate factorization (AF) 
algorithm for solution of the unsteady transonic small-disturbance equation including entropy 
and vorticity effects. The AF algorithm has been shown to be very efficient for steady or 
unsteady transonic flow problems. It can provide accurate solutions in only several hundred 
time steps yielding a significant computational cost savings when compared to alternative 
methods. For reasons of practicality and affordability, an efficient algorithm and a fast 
computer code are requirements for realistic aircraft applications. 

Results were presented for several configurations including the General Dynamics one-ninth 
scale F-16C aircraft model and the ONERA M6 wing, to demonstrate various capabilities of the 
CAP-TSD code. For the F-16C aircraft model, calculated steady pressure distributions 
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compared well with the experimental data. Unsteady pressures for the entire F-16C aircraft 
undergoing a rigid pitching motion were also presented. Comparisons with parallel wing alone 
results revealed aerodynamic interference effects of the additional aircraft components on wing 
unsteady pressures. These effects emphasize the importance of including all components in the 
calculation. The CAP-TSD code thus provides the capability of modeling complete aircraft 
configurations for realistic transonic unsteady aerodynamic and aeroelastic analyses. For the 
ONERA M6 wing, CAP-TSD results were presented both with and without entropy and vorticity 
effects. The results obtained by including these effects were in very good agreement with Euler 
results in terms of predicting the shock location and strength. The CAP-TSD code includes the 
entropy and vorticity effects while retaining the relative simplicity and cost efficiency of the 
small-disturbance potential formulation. Therefore, the capability provides the aeroelastician 
with an affordable method to analyze relatively difficult transonic cases. 

Results were also presented from a flutter analysis of a 45° sweptback wing. The flutter 
boundaries from CAP-TSD (linear) were in agreement with parallel subsonic linear theory 
results and compared well with the experimental flutter data for subsonic and supersonic 
freestream Mach numbers. The nonlinear CAP-TSD flutter results also compared favorably 
with the experimental data which is the first step toward validating the code for general 
transonic aeroelastic analysis of more complex configurations. 
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CAP-TSD modeling of the General Dynamics one-ninth scale F-16C aircraft model. 
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Comparison between CAP-TSD and experimental steady pressure distributions on the 
wing and tail of the F-16C aircraft model at M = 0.9 and ao = 2.38°. 
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CAP-TSD unsteady pressure distributions on the wing and tail upper surfaces of the 
F-16C aircraft model due to complete airplane rigid pitching at M = 0.9, ao = 2.38°, 
ai = 0.5°, and k = 0.1. 
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Comparison of steady pressure distributions for the ONERA M6 wing at M = 0.92 and 
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Fig. 5 Planview of 45 ° sweptback wing. 
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Comparisons between calculated and experimental flutter speed index of the 45 
sweptback wing. 
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